The previous two questions are similar. Fitting the models in JAGS is not that challenging. The interesting parts of these problems were selecting and justifying the prior and interpreting the posterior.

There is some lattitude regarding how to select a prior. Let’s suppose the likelihood is specified by

\[Y_{i|g} \sim {\sf Norm}(\mu_g, \sigma_g)\]

Then we need priors for \(\mu_g\) and \(\sigma_g\). (Technically, we need a joint prior for \(\mu_1\), \(\mu_2\), \(\sigma_1\) and \(\sigma_2\). But if the marginal priors are independent, then specifying the marginal priors determines the joint prior.)

Some important things about selecting the prior:

• A normal prior is a fine choice for \(\mu_g\).

• The mean of this normal prior for \(\mu_g\) should be near where we believe the mean height of people (in one group) is (prior to using our data).
  0 is not a good choice in this situation unless you use a very large \(\sigma\) because we don’t expect the mean (for eiether group) to be close to 0.

• It is fine to use the same mean for the priors for both genders. The data will push those apart in the posterior as long as the prior isn’t too narrow.

• The standard deviation/precision of the prior for \(\mu_g\) is not a measure of variability in the population. It is a measure of your certainty about the mean you are providing to the prior. Your explanation of how you chose the prior for the standard deviation should reflect that you understand this distinction.

• You also need a prior for the standard deviation \(\sigma_g\). \(\sigma_g\) represents the amount of variability of \(y\) (height or pulse) in each gender group. This is the place to encorporate what you know about how large that standard deviation might be.

Note the use of separate standard deviations is not an “assumption that they are different”, but it allows them to be different and so gives us a way to test wether they are the same or different (as you were asked to do in part c).

Don’t forget to do some diagnostics. Some of you made small errors/typos that led to models that were clearly not correct but didn’t notice.

A note on the posterior: There is only one posterior (and only one prior, and one likelihood). Typically the posterior is a joint distribution, so we can talk about the marginal distributions of its components. We will often refer to these as the posterior for \(\mu\) or the posterior for \(\sigma\). It would be more explicit to refer to these as the marginal posteriors for these quantities.

When checking whether two parameters are different, look at the posterior distribution for the difference, not the amount of overlap in the two posterior distributions.

If the posterior distribution of a difference straddles 0, this does not allow you to conclude that there is no difference. If you use a ROPE and the entire HDI is contained in the ROPE, then you can conclude there is no practical difference (at least to the degree specified by the probability of your HDI).

The two groups look different, but neither one looks particularly symmetric. (Other types of plots could be used here as well.)

gf_violin(DaysLive ~ Group, data = RatLives)

Nevertheless, we will proceed with this model for the moment. Let’s build a generic version of the model in JAGS that we can recycle for other purposes as we go along.

# y ~ g where y is metric and g is 1 or 2

metric_dich_model <- function() {
  for (i in 1:Nobs) {
    y[i] ~ dt(mu[g[i]], 1 / sigma[g[i]]^2, nu[g[i]])
  }
  for (g in 1:2) {
    mu[g]         ~ dnorm(mu_mean, 1 / mu_sd^2)
    sigma[g]      ~ dunif(sigma_lo, sigma_hi)
    nuMinusOne[g] ~ dexp(1 / 29.0)
    nu[g]        <- nuMinusOne[g] + 1
    log10nu[g]   <- log(nu[g]) / log(10)
  }
  delta_mu     <- mu[2] - mu[1]
  delta_sqrtmu <- sqrt(abs(mu[2])) - sqrt(abs(mu[1]))  # for sake of second attempt
  sigma_ratio  <- sigma[2] / sigma[1]
}

Now let’s fit the model.

rats_jags <-   
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = nrow(RatLives),
      y = RatLives$DaysLive,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700,      # approximate average life of a rat
      mu_sd    = 300,      # pretty sure mean is between 100 and 1300 (for each group)
      sigma_lo = 60 / 1000, 
      sigma_hi = 60 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 4,
    n.thin = 2,
    DIC = FALSE
  )

Let’s take a look a what we have.

rats_jags

## Inference for Bugs model at "/var/folders/py/txwd26jx5rq83f4nn0f5fmmm0000gn/T//RtmpqhObxc/model1610a59ae5540.txt", fit using jags,
##  4 chains, each with 5000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 8000 iterations saved
##              mu.vect sd.vect    2.5%     25%      50%      75%    97.5%
## delta_mu     298.101  33.927 229.138 275.299  299.174  321.839  360.356
## log10nu[1]     0.484   0.163   0.199   0.373    0.469    0.581    0.842
## log10nu[2]     0.874   0.404   0.335   0.573    0.771    1.105    1.834
## mu[1]        706.029  10.643 684.321 698.888  706.240  713.307  726.562
## mu[2]       1004.129  32.369 937.931 982.560 1005.769 1026.805 1063.531
## sigma[1]      79.650  11.628  58.556  71.476   79.118   87.317  104.168
## sigma[2]     232.669  36.021 167.582 205.475  231.624  259.367  302.287
## sigma_ratio    2.986   0.653   1.879   2.525    2.925    3.387    4.423
##              Rhat n.eff
## delta_mu    1.001  8000
## log10nu[1]  1.001  5300
## log10nu[2]  1.002  2500
## mu[1]       1.001  8000
## mu[2]       1.001  8000
## sigma[1]    1.001  5600
## sigma[2]    1.001  8000
## sigma_ratio 1.001  8000
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

rats_mcmc <- as.mcmc(rats_jags)
mcmc_combo(rats_mcmc, regex_pars = c("delta", "ratio", "log"))
diag_mcmc(rats_mcmc)

Convergence looks good, but \(\sigma\) and \(\nu\) are being inefficiently sampled. (If we were really going to use the model, we would like run longer to boost the effective sample size.)

Suppse we are only interested in differences in average lifetime of at least 60 days, and differences in \(\sigma\) of at least 10%.

plot_post(posterior(rats_jags)$delta_mu, ROPE = c(-60, 60), xlab = "delta_mu")

## $posterior
##           ESS     mean   median     mode
## var1 2070.804 298.1006 299.1735 299.4048
## 
## $hdi
##   prob     lo       hi
## 1 0.95 232.52 362.6606
## 
## $ROPE
##    lo hi P(< ROPE) P(in ROPE) P(> ROPE)
## 1 -60 60         0          0         1

hdi(posterior(rats_jags), regex_pars = "delta_mu")

##        par     lo       hi prob
## 1 delta_mu 232.52 362.6606 0.95

plot_post(posterior(rats_jags)$sigma_ratio, ROPE = c(0.9, 1.1), xlab = "sigma_ratio")

## $posterior
##           ESS     mean   median     mode
## var1 1831.936 2.986193 2.924575 2.731553
## 
## $hdi
##   prob     lo       hi
## 1 0.95 1.7616 4.263581
## 
## $ROPE
##    lo  hi P(< ROPE) P(in ROPE) P(> ROPE)
## 1 0.9 1.1         0          0         1

hdi(posterior(rats_jags), regex_pars = "sigma_ratio")

##           par     lo       hi prob
## 1 sigma_ratio 1.7616 4.263581 0.95

hdi(posterior(rats_jags), regex_pars = "log")

##         par        lo        hi prob
## 1 log10nu.1 0.1813248 0.8099559 0.95
## 2 log10nu.2 0.2839673 1.7514411 0.95

According to this model,

• There is a clear difference in mean longevity, with all of the posterior distribution for the difference in means is well outside a ROPE of \(\pm 2\) months (60 days).
  Indeed, the rats on the restricted diet seem to live 8 months to a year longer on average.
• the standard deviations also appear to be different.
• group 1 (ad lib) seems to have clearly heavier than normal tails.
• group 2 (restricted) may also have heavier tails but the evidence is not as clear (much of the posterior credibility of \(log_10(\nu)\) is less than 1.5).

But modeling the distribution of the response as a symmetric (t) distribution doesn’t seem to align with the data. So let’s try squaring the response.

gf_violin((DaysLive^2) ~ Group, data = RatLives) 
gf_dhistogram(~ (DaysLive^2) | Group ~ ., data = RatLives, bins = 40) %>% 
  gf_fitdistr()

The squared lives do appear to be more symmetrically distributed.

Let’s a fit a new model wit the lifetimes squared. Note that the mean of a squared distribuion is not the mean of the original distribution squared:

\[E(X^2) = E(X)^2 + Var(X) \neq E(X)^2\]

rats2_jags <-   
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "delta_sqrtmu", "log10nu"),
    data = list(
      Nobs = nrow(RatLives),
      y = RatLives$DaysLive^2,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700^2 + 300^2,
      mu_sd    = (700^2 + 300^2) / 3,  # takes 3 sd's to get to negative values
      sigma_lo = 500 / 1000, 
      sigma_hi = 500 * 1000
    ),
    n.iter   = 5000,
    n.burnin = 1000,
    n.chains = 4,
    n.thin   = 2,
    DIC      = FALSE
  )

summary(rats2_jags)

## fit using jags
##  4 chains, each with 5000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 8000 iterations saved
##                  mu.vect   sd.vect       2.5%        25%         50%
## delta_mu      508736.608 48178.970 415019.324 476095.302  508838.348
## delta_sqrtmu     298.580    25.430    248.680    281.440     298.647
## log10nu[1]         1.030     0.391      0.431      0.727       0.969
## log10nu[2]         1.530     0.310      0.911      1.313       1.541
## mu[1]         492920.082 15987.963 460788.700 482564.279  492914.583
## mu[2]        1001656.690 45406.261 911951.828 971722.640 1001280.963
## sigma[1]      132907.193 17554.461  99642.494 120919.095  133018.534
## sigma[2]      453425.474 26875.799 396808.997 435422.322  455529.199
## sigma_ratio        3.474     0.526      2.621      3.098       3.412
##                      75%       97.5%  Rhat n.eff
## delta_mu      540871.192  603009.022 1.002  2700
## delta_sqrtmu     315.598     347.809 1.002  2800
## log10nu[1]         1.293       1.876 1.001  8000
## log10nu[2]         1.759       2.085 1.002  2700
## mu[1]         503630.988  523900.152 1.002  2500
## mu[2]        1032322.489 1090053.088 1.002  2400
## sigma[1]      144970.909  167127.684 1.001  8000
## sigma[2]      473956.413  496266.281 1.002  2500
## sigma_ratio        3.787       4.649 1.001  8000

rats2_mcmc <- as.mcmc(rats2_jags)
plot_post(posterior(rats2_jags)$sigma_ratio, xlab = "sigma_ratio")

## $posterior
##           ESS     mean   median     mode
## var1 2237.822 3.474205 3.411965 3.345158
## 
## $hdi
##   prob       lo       hi
## 1 0.95 2.473909 4.463433

plot_post(posterior(rats2_jags)$delta_sqrtmu, xlab = "delta_sqrtmu")

## $posterior
##           ESS     mean   median   mode
## var1 6704.354 298.5803 298.6468 299.72
## 
## $hdi
##   prob       lo       hi
## 1 0.95 249.4157 348.2196

In this model we have strong evidence that the standard deviations (on the squared scale) and the means are different. The distributions of squared lifetimes appear much more normal than the untransformed lifetimes.

Let’s compare the HDIs for the difference in mean lifetimes for the two models:

HDIs <- 
  bind_rows(
    hdi(posterior(rats_jags), regex_pars = "delta_mu") %>% 
      mutate(model = "model 1"),
    hdi(posterior(rats2_jags), regex_pars = "delta_sqrtmu") %>% 
      mutate(model = "model 2")
  ) %>% 
  select(model, matches("*"))  # get model to the first column
HDIs

##     model          par       lo       hi prob
## 1 model 1     delta_mu 232.5200 362.6606 0.95
## 2 model 2 delta_sqrtmu 249.4157 348.2196 0.95

gf_dens(~delta_mu, data = posterior(rats_jags), color = ~ "model 1") %>%
  gf_dens(~delta_sqrtmu, data = posterior(rats2_jags), color = ~ "model 2") %>%
  gf_vline(xintercept = ~lo, color = ~model, data = HDIs) %>%
  gf_vline(xintercept = ~hi, color = ~model, data = HDIs) %>%
  gf_labs(x = "difference in means", color = "model") 

The centers of the posterior distributions are quite similar, but model 2 gives us a narrower HDI for the difference in the means.

rats_prior1 <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = 0,
      y = RatLives$DaysLive,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700,
      mu_sd    = 300,
      sigma_lo = 60 / 1000, 
      sigma_hi = 60 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )
rats_prior2 <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "delta_sqrtmu", "log10nu"),
    data = list(
      Nobs = 0,
      y = RatLives$DaysLive^2,
      g = as.numeric(factor(RatLives$Group)),
      mu_mean  = 700^2 + 300^2,
      mu_sd    = (700^2 + 300^2) / 3,         
      sigma_lo = 500 / 1000, 
      sigma_hi = 500 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )

gf_dens(~ delta_mu, data = posterior(rats_prior1), color = ~"model 1") %>%
gf_dens(~ delta_sqrtmu, data = posterior(rats_prior2), color = ~"model 2")

The prior for the difference in the mean lifetimes is a bit more spread out in the first model and less for the second model. If our goals was keep them similar, we should adjust one or the other. Despite this difference, both models agree that delta_mu is approximately 300. (Both priors consider 300 to be a plausible value, so it isn’t surprising that they agree in the posterior.)

flies_jags <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = nrow(ShohatOphirKAMH2012dataReduced),
      y = ShohatOphirKAMH2012dataReduced$PreferenceIndex,
      g = as.numeric(factor(ShohatOphirKAMH2012dataReduced$Group)),
      mu_mean  = 0,      # 0 = no preference
      mu_sd    = 0.5,    # range is (-1, 1), this will sample all of that range, but prefer near 0
      sigma_lo = 0.0001, # close to 0
      sigma_hi = 1       # 1 is max possible because values in (-1, 1) 
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )

## module glm loaded

## module dic loaded

(Omitting diagnositics to save space.)

Posterior distribution and HDI for delta_mu and effect size.

flies_post <- 
  posterior(flies_jags) %>% 
  mutate(efsize = delta_mu / sqrt((sigma.1^2 + sigma.2^2) / 2))
gf_dens( ~ efsize, data = flies_post)
gf_dens( ~ delta_mu, data = flies_post)
hdi(flies_post, pars = c("efsize", "delta_mu"))

##        par         lo        hi prob
## 1 delta_mu 0.08399878 0.4945389 0.95
## 2   efsize 0.25500843 2.3855492 0.95

Conclusion:

• There appears to be a credible difference in preference (we would reject the hypothesis that there is no difference).
• The effect size is credibly bigger than 0.2 (a benchmark often used in pyschology to indicate a “big” effect).

flies2_jags <-
  jags(
    model = metric_dich_model,
    parameters.to.save = 
      c("mu", "sigma", "delta_mu", "sigma_ratio", "log10nu"),
    data = list(
      Nobs = nrow(ShohatOphirKAMH2012dataReduced),
      y = ShohatOphirKAMH2012dataReduced$GrandTotal,
      g = as.numeric(factor(ShohatOphirKAMH2012dataReduced$Group)),
      mu_mean  = 150,    # rough order of maginitude of amonts eaten
      mu_sd    = 100,    # very broad
      sigma_lo = 10 / 1000,
      sigma_hi = 10 * 1000
    ),
    n.iter = 5000,
    n.burnin = 1000,
    n.chains = 2,
    DIC = FALSE
  )

Posterior distribution and HDI for delta_mu and effect size.

flies2_post <- 
  posterior(flies2_jags) %>% 
  mutate(efsize = delta_mu / sqrt((sigma.1^2 + sigma.2^2) / 2))
gf_dens( ~ efsize, data = flies2_post)
gf_dens( ~ delta_mu, data = flies2_post)
hdi(flies2_post, pars = c("efsize", "delta_mu"))

##        par          lo        hi prob
## 1 delta_mu -19.2292020 22.820278 0.95
## 2   efsize  -0.8284769  0.991917 0.95

Conclusion:

• There is a lot of variability, but neither group clearly eats more than the other.
• Both large and small effect sizes are credible.
• We cannot attribute the preference difference in the previous part to differences in total consumption. (If there were a clear indication that one group consumed more, that could explain why the preference ratio was different from 0.)

set.seed(54321)
library(dplyr)
GaltonM <-
  mosaicData::Galton %>% 
  filter(sex == "M") %>%
  group_by(family) %>%
  sample_n(1) %>%
  mutate(midparent = (father + mother) / 2)

galton_model <- function() {
  for (i in 1:length(y)) {
    y[i]   ~ dt(mu[i], 1/sigma^2, nu)
    mu[i] <- alpha0 + alpha1 * (x[i] - mean(x))
  }
  sigma ~ dunif(6/100, 6 * 100)
  nuMinusOne ~ dexp(1/29)
  nu <- nuMinusOne + 1
  alpha0 ~ dnorm(0, 1/100^2)   # 100 is order of magnitude of data
  alpha1 ~ dnorm(0, 1/4^2)     # expect roughly 1-1 slope
  beta1 <- alpha1
  beta0 <- alpha0 - mean(x) * alpha1
}

library(R2jags)
library(mosaic)
galton_jags <-
  jags(
    model = galton_model,
    data = list(y = GaltonM$height, x = GaltonM$midparent),
    parameters.to.save = c("beta0", "beta1", "sigma", "nu"),
    n.iter   = 5000,
    n.burnin = 2000,
    n.chains = 4,
    n.thin   = 1
  )

Post <- 
  posterior(galton_jags) %>% 
  mutate(
    avg_ht = beta0 + beta1 * 68.5,
    ind_ht = avg_ht + rstudent_t(12000, mean = 0, sigma = sigma, df = nu))

Post %>% hdi(0.95, pars = "avg_ht")

##      par       lo       hi prob
## 1 avg_ht 70.17743 71.15686 0.95

Post %>% hdi(0.50, pars = "ind_ht")

##      par       lo       hi prob
## 1 ind_ht 69.20809 72.18985  0.5

gf_dens(~avg_ht, data = Post, color = ~ "Average") %>% 
  gf_dens(~ind_ht, data = Post, color = ~ "Individual") %>%
  gf_labs(x = "height of male with midparent = 68.5")
 
# plot_post(Post$avg_ht, hdi_prob = 0.95, xlab = "average")
# plot_post(Post$ind_ht, hdi_prob = 0.50, xlab = "individual")

Centering on x was needed to break the correlation between the slope and intercept (which makes JAGS much more efficient). It can sometimes make the intercept easier to interpret as well.

Centering on y does not hurt anything, but the correlation has already been removed by centering on x, so centering y isn’t necessary. If we do center y, then our response essentially becomes how different is y from the average y in our data.

Note: We don’t need to center on the mean. Sometimes centering is done with convenient benchmark value. As long as the benchmark is near the mean of the values in our data, correlation will be greatly reduced. Said another way: we should pick a benchmark that is imidst the data – or design the study to collect data centered roughly on the benchmark value.

Let’s begin by fitting the models. (Be sure to use results = "hide" in your R chunk header to avoid including all of the progress indicators in your output. You may also like to use cache = TRUE so that the model isn’t refit each time you knit the document.)

model1_brm <-
  brm(
    bf(longevity ~ group + thorax, sigma ~ group + thorax),
    data = FruitflyReduced, 
    family = student(),
  )

## Compiling the C++ model

## Start sampling

model2_brm <-
  brm(
    bf(longevity ~ group, sigma ~ group),
    data = FruitflyReduced, 
    family = student(),
  )

## Compiling the C++ model

## recompiling to avoid crashing R session

## Start sampling

post1 <- posterior(model1_brm)
post2 <- posterior(model2_brm)

a The posterior for b_thorax in model 1 is well separated from 0, and the 95% HDI is

    h <- hdi(post1, pars = "b_thorax"); h

##        par       lo       hi prob
## 1 b_thorax 117.3338 158.1048 0.95

The units here is days/mm. So roughly for every 0.1 mm longer their thorax is, flies live on average between 11.7 and 15.8 days longer. (For example, the average longevity for flies with 0.8 mm thorax is roughly 11.7 to 15.8 days longer than for flies with 0.7mm).

b We can construct a contrast to measure the difference of interest. For model 1 this looks like

\[\begin{align*} \frac{\mu_{V1} + \mu_{V8}}{2} - \frac{\mu_{P1} + \mu_{P8} + \mu_{N}}{3} &= \frac{\beta_{V1} + \beta_0 + \beta_t t + \beta_{V8} + \beta_0 + \beta_t t}{2} - \frac{\beta_{P1} + \beta_0 + \beta_t t + \beta_{P8} + \beta_0 + \beta_t t + \beta_0 + \beta_t t}{3} \\ &= \frac{\beta_{V1} + \beta_{V8}}{2} - \frac{\beta_{P1} + \beta_{P8}}{3} \end{align*}\]

This will be the same for model 2 since it doesn’t have the thorax terms to begin with.

    post1 <- post1 %>%
      mutate(VvsOther = 
               (b_groupVirgin1 + b_groupVirgin8) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    post2 <- post2 %>%
      mutate(VvsOther = 
               (b_groupVirgin1 + b_groupVirgin8) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    bind_rows(
      hdi(post1, pars = "VvsOther") %>% mutate(model = 1),
      hdi(post2, pars = "VvsOther") %>% mutate(model = 2)
    )

##        par        lo        hi prob model
## 1 VvsOther -18.82952 -11.51699 0.95     1
## 2 VvsOther -21.49390 -10.71867 0.95     2

c

There is a lot of overlap of the two HDIs, but the HDI from model 1 is narrower (and shifted a bit lower). So model 1 provides strong evidence that the mean lifetime for the two groups differ (on average).

    post1 <- post1 %>%
      mutate(V1vsnoV = 
               (b_groupVirgin1) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    post2 <- post2 %>%
      mutate(V1vsnoV = 
               (b_groupVirgin1) / 2 -
               (b_groupPregnant1 + b_groupPregnant8) / 3
      )
    bind_rows(
      hdi(post1, pars = "V1vsnoV") %>% mutate(model = 1),
      hdi(post2, pars = "V1vsnoV") %>% mutate(model = 2)
    )

##       par        lo         hi prob model
## 1 V1vsnoV -8.835859 -2.5139140 0.95     1
## 2 V1vsnoV -8.026582  0.9415611 0.95     2

d

    post1 <- post1 %>% 
      mutate(
        SigmaP1T7 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept + b_sigma_thorax * 0.7),
        SigmaP1T8 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept + b_sigma_thorax * 0.8),
        SigmaP1T9 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept + b_sigma_thorax * 0.9),
        SigmaV8T7 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept + b_sigma_thorax * 0.7),
        SigmaV8T8 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept + b_sigma_thorax * 0.8),
        SigmaV8T9 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept + b_sigma_thorax * 0.9)
        )
    post2 <- post2 %>% 
      mutate(
        SigmaP1T7 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept),
        SigmaP1T8 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept),
        SigmaP1T9 = exp(b_sigma_groupPregnant1 + b_sigma_Intercept),
        SigmaV8T7 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept),
        SigmaV8T8 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept),
        SigmaV8T9 = exp(b_sigma_groupVirgin8 + b_sigma_Intercept)
        )
    Sigmas <-
      bind_rows(
        hdi(post1, pars = grep("Sigma", names(post1), value = TRUE)) %>% mutate(model = 1),
        hdi(post2, pars = grep("Sigma", names(post2), value = TRUE)) %>% mutate(model = 2)
      ) %>% 
      arrange(par)
    Sigmas

##          par        lo        hi prob model
## 1  SigmaP1T7  6.264782 14.171235 0.95     1
## 2  SigmaP1T7 11.243253 20.607381 0.95     2
## 3  SigmaP1T8  8.486591 15.818850 0.95     1
## 4  SigmaP1T8 11.243253 20.607381 0.95     2
## 5  SigmaP1T9  9.589138 19.396911 0.95     1
## 6  SigmaP1T9 11.243253 20.607381 0.95     2
## 7  SigmaV8T7  3.748734  8.198663 0.95     1
## 8  SigmaV8T7  8.794834 16.404818 0.95     2
## 9  SigmaV8T8  4.999383  9.217360 0.95     1
## 10 SigmaV8T8  8.794834 16.404818 0.95     2
## 11 SigmaV8T9  5.772719 11.641738 0.95     1
## 12 SigmaV8T9  8.794834 16.404818 0.95     2

    Sigmas %>%
      gf_errorbar(lo + hi ~ par, color = ~factor(model), 
                  position = "dodge") %>%
      gf_labs(color = "model", title = "85% HPDIs for sigma")

Since thorax size explains some of the variability in lifetime, the standard deviation estimates of model 1 are smaller than those for model 2. The width of the HDI also appears smaller.

e Inclduing thorax seems to be a good idea. Longevity is positively associated with thorax size, and including thorax in the model adjusts for any possible differences in the sizes of the flies in the different groups. The smaller estimates of \(\sigma\) indicate that the model with thorax is able to give for precises estimates.

a Looking at the data:

b Fitting a model

sw_brm <- brm(SeaweedAmt ~ Zone * Grazer, data = Seaweed)

## Compiling the C++ model

## Start sampling

Since we only have two observations in each zone/grazer group, it would be very difficult to estimate variances in each group. Since the interation model has a mean for each zone/grazer combo, if we added a separate variance for each the model would fit the data too well: If there are only two values and we know the mean and variance, we can recover the two values.

c

The biggest hassle here is creating the expression for the contrast you are interested in. I’ve written a little function to simplify this a bit, but I’ve hidden the code.

Little fish vs no little fish (no big fish in either case)

c1 <- 
  create_contrast(M[,"None"], M[, "f"], 
                  a_denom = 8, b_denom = 8, rel = "=")
c2 <- 
  create_contrast(M[,"L"], M[, "Lf"], 
                  a_denom = 8, b_denom = 8, rel = "=")
c3 <- 
  create_contrast(M[,c("L", "None")], M[, c("Lf", "f")], 
                  a_denom = 16, b_denom = 16, rel = "=")
c1

## [1] "( GrazerNone + ZoneB + GrazerNone + ZoneB:GrazerNone + ZoneC + GrazerNone + ZoneC:GrazerNone + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneE + GrazerNone + ZoneE:GrazerNone + ZoneF + GrazerNone + ZoneF:GrazerNone + ZoneG + GrazerNone + ZoneG:GrazerNone + ZoneH + GrazerNone + ZoneH:GrazerNone ) / 8 = (  + ZoneB + ZoneC + ZoneD + ZoneE + ZoneF + ZoneG + ZoneH ) / 8"

c2

## [1] "( GrazerL + ZoneB + GrazerL + ZoneB:GrazerL + ZoneC + GrazerL + ZoneC:GrazerL + ZoneD + GrazerL + ZoneD:GrazerL + ZoneE + GrazerL + ZoneE:GrazerL + ZoneF + GrazerL + ZoneF:GrazerL + ZoneG + GrazerL + ZoneG:GrazerL + ZoneH + GrazerL + ZoneH:GrazerL ) / 8 = ( GrazerLf + ZoneB + GrazerLf + ZoneB:GrazerLf + ZoneC + GrazerLf + ZoneC:GrazerLf + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneE + GrazerLf + ZoneE:GrazerLf + ZoneF + GrazerLf + ZoneF:GrazerLf + ZoneG + GrazerLf + ZoneG:GrazerLf + ZoneH + GrazerLf + ZoneH:GrazerLf ) / 8"

c3

## [1] "( GrazerL + ZoneB + GrazerL + ZoneB:GrazerL + ZoneC + GrazerL + ZoneC:GrazerL + ZoneD + GrazerL + ZoneD:GrazerL + ZoneE + GrazerL + ZoneE:GrazerL + ZoneF + GrazerL + ZoneF:GrazerL + ZoneG + GrazerL + ZoneG:GrazerL + ZoneH + GrazerL + ZoneH:GrazerL + GrazerNone + ZoneB + GrazerNone + ZoneB:GrazerNone + ZoneC + GrazerNone + ZoneC:GrazerNone + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneE + GrazerNone + ZoneE:GrazerNone + ZoneF + GrazerNone + ZoneF:GrazerNone + ZoneG + GrazerNone + ZoneG:GrazerNone + ZoneH + GrazerNone + ZoneH:GrazerNone ) / 16 = ( GrazerLf + ZoneB + GrazerLf + ZoneB:GrazerLf + ZoneC + GrazerLf + ZoneC:GrazerLf + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneE + GrazerLf + ZoneE:GrazerLf + ZoneF + GrazerLf + ZoneF:GrazerLf + ZoneG + GrazerLf + ZoneG:GrazerLf + ZoneH + GrazerLf + ZoneH:GrazerLf +  + ZoneB + ZoneC + ZoneD + ZoneE + ZoneF + ZoneG + ZoneH ) / 16"

hypothesis(sw_brm, c1) 

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerNone+Zone... = 0     9.18      3.47     2.47    16.18         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(sw_brm, c2) 

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerL+ZoneB+G... = 0     2.73      3.38    -4.07      9.4         NA
##   Post.Prob Star
## 1        NA     
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(sw_brm, c3) 

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerL+ZoneB+G... = 0     5.95      2.46     1.16    10.73         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(sw_brm, c1) %>% plot()
hypothesis(sw_brm, c2) %>% plot()
hypothesis(sw_brm, c3) %>% plot()

Limpets vs No Limpets

c4 <- 
  create_contrast(
    M[, c("None", "f", "fF")], M[, c("L", "Lf", "LfF")], 
    a_denom = 24, b_denom = 24, rel = "=")
c4

## [1] "( GrazerNone + ZoneB + GrazerNone + ZoneB:GrazerNone + ZoneC + GrazerNone + ZoneC:GrazerNone + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneE + GrazerNone + ZoneE:GrazerNone + ZoneF + GrazerNone + ZoneF:GrazerNone + ZoneG + GrazerNone + ZoneG:GrazerNone + ZoneH + GrazerNone + ZoneH:GrazerNone +  + ZoneB + ZoneC + ZoneD + ZoneE + ZoneF + ZoneG + ZoneH + GrazerfF + ZoneB + GrazerfF + ZoneB:GrazerfF + ZoneC + GrazerfF + ZoneC:GrazerfF + ZoneD + GrazerfF + ZoneD:GrazerfF + ZoneE + GrazerfF + ZoneE:GrazerfF + ZoneF + GrazerfF + ZoneF:GrazerfF + ZoneG + GrazerfF + ZoneG:GrazerfF + ZoneH + GrazerfF + ZoneH:GrazerfF ) / 24 = ( GrazerL + ZoneB + GrazerL + ZoneB:GrazerL + ZoneC + GrazerL + ZoneC:GrazerL + ZoneD + GrazerL + ZoneD:GrazerL + ZoneE + GrazerL + ZoneE:GrazerL + ZoneF + GrazerL + ZoneF:GrazerL + ZoneG + GrazerL + ZoneG:GrazerL + ZoneH + GrazerL + ZoneH:GrazerL + GrazerLf + ZoneB + GrazerLf + ZoneB:GrazerLf + ZoneC + GrazerLf + ZoneC:GrazerLf + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneE + GrazerLf + ZoneE:GrazerLf + ZoneF + GrazerLf + ZoneF:GrazerLf + ZoneG + GrazerLf + ZoneG:GrazerLf + ZoneH + GrazerLf + ZoneH:GrazerLf + GrazerLfF + ZoneB + GrazerLfF + ZoneB:GrazerLfF + ZoneC + GrazerLfF + ZoneC:GrazerLfF + ZoneD + GrazerLfF + ZoneD:GrazerLfF + ZoneE + GrazerLfF + ZoneE:GrazerLfF + ZoneF + GrazerLfF + ZoneF:GrazerLfF + ZoneG + GrazerLfF + ZoneG:GrazerLfF + ZoneH + GrazerLfF + ZoneH:GrazerLfF ) / 24"

hypothesis(sw_brm, c4) 

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerNone+Zone... = 0    28.22      1.98    24.29    32.12         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(sw_brm, c4) %>% plot()

Zone A vs Zone D

c5 <- 
  create_contrast( M["A", ], M["D", ],
    a_denom = 6, b_denom = 6, rel = "=")
c5

## [1] "(  + GrazerfF + GrazerL + GrazerLf + GrazerLfF + GrazerNone ) / 6 = ( ZoneD + ZoneD + GrazerfF + ZoneD:GrazerfF + ZoneD + GrazerL + ZoneD:GrazerL + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneD + GrazerLfF + ZoneD:GrazerLfF + ZoneD + GrazerNone + ZoneD:GrazerNone ) / 6"

hypothesis(sw_brm, c5) 

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((+GrazerfF+Graze... = 0   -45.53      3.97   -53.22   -37.77         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(sw_brm, c5) %>% plot()

c6 <- 
  create_contrast( 
    c(M["A", c("L", "Lf", "LfF")], M["D", c("None", "f", "fF")]),
    c(M["D", c("L", "Lf", "LfF")], M["A", c("None", "f", "fF")]),
    a_denom = 6, b_denom = 6, rel = "=")
c6

## [1] "( GrazerL + GrazerLf + GrazerLfF + ZoneD + GrazerNone + ZoneD:GrazerNone + ZoneD + ZoneD + GrazerfF + ZoneD:GrazerfF ) / 6 = ( ZoneD + GrazerL + ZoneD:GrazerL + ZoneD + GrazerLf + ZoneD:GrazerLf + ZoneD + GrazerLfF + ZoneD:GrazerLfF + GrazerNone +  + GrazerfF ) / 6"

hypothesis(sw_brm, c6) 

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 ((GrazerL+GrazerL... = 0    20.55      3.96    12.64    28.33         NA
##   Post.Prob Star
## 1        NA    *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(sw_brm, c6) %>% plot()

model1_brm <- brm(Y ~ Group, data = NonhomogVar)

## Compiling the C++ model

## Start sampling

model2_brm <- brm(bf(Y ~ Group, sigma ~ Group), data = NonhomogVar)

## Compiling the C++ model
## Start sampling

The second of these models uses a different value of \(\sigma\) for each group. The first uses the same \(\sigma\) for all the groups.

Priors:

• intercept (both models): \({\sf T}(3, 100, 10)\)

  df_stats(~Y, data = NonhomogVar, mean, sd)

  ##   mean_Y     sd_Y
  ## 1  100.5 6.228965

  So this is using a (order-of-magnitude rounded) mean and standard deviation of the response to determine the prior.

• \(\log(\sigma)\) (model 1): \({\sf T}(3, 0, 10)\)
• \(\log(\sigma_A)\) (model 2): \({\sf T}(3, 0, 10)\)

• other linear model coefficients (including the offsets for \(\log(\sigma)\) in model 2): improper uniform prior

As expected, the model that fits separate \(\sigma\)’s for each group estimates two smaller values and to larger values. Model 1 has a single \(\sigma\) that falls between the two.

Post1 <- posterior(model1_brm)
Post2 <- posterior(model2_brm) %>%
  mutate(
    sigmaA = exp(b_sigma_Intercept),     # don't forget exp()!
    sigmaB = exp(b_sigma_Intercept + b_sigma_GroupB),
    sigmaC = exp(b_sigma_Intercept + b_sigma_GroupC),
    sigmaD = exp(b_sigma_Intercept + b_sigma_GroupD)
  ) 
Post2a <- 
  Post2 %>% 
  select(sigmaA, sigmaB, sigmaC, sigmaD) %>%
  tidyr::gather(group, sigma, sigmaA:sigmaD)

gf_dens(~sigma, data = Post1, size = 0.7) %>%
  gf_dens(~sigma, data = Post2a, color = ~group, size = 1, alpha = 0.5)

Also as expected, Model 2 finds no credible difference between the standard deviations of group B and group C, but is quite convinced that there is a difference bewteen groups A and B and between (A or D) vs (B or C).

gf_dens( ~(sigmaA - sigmaB), data = Post2, color = ~"A vs B", size = 1) %>%
gf_dens( ~(sigmaB - sigmaC), data = Post2, color = ~"B vs C", size = 1) %>%
gf_dens(~((sigmaA + sigmaD) / 2 - (sigmaB + sigmaC) / 2), color = ~ "AD vs BC",
        data = Post2, size = 1)

The HDIs tell the same story as the plots above.

bind_rows(
  hdi(Post2, pars =  "sigma: A vs B" ~ (sigmaA - sigmaB)),
  hdi(Post2, pars =  "sigma: B vs C" ~ (sigmaB - sigmaC)),
  hdi(Post2, 
      pars =  "sigma: AD vs BC" ~ ((sigmaA + sigmaD) / 2 - (sigmaB + sigmaC) / 2))
)

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

##                 par         lo       hi prob
## 1   "sigma: A vs B"  4.9898263 9.996451 0.95
## 2   "sigma: B vs C" -0.4389844 0.434171 0.95
## 3 "sigma: AD vs BC"  5.6391593 9.116402 0.95

Looking at how the two models compare the means we see that the results are quite similar when we compare group A to group B, but when we compare group B to group C, Model 2 finds compelling evidence that \(\mu_B < \mu_C\) while Model 1 does not. The reason for the difference is the Model 1 estimates each group to have a much larger standard deviation (influenced by groups A and D). The extra modeled variability within the groups hides the difference between the two groups.

bind_rows(
  hdi(Post1, pars =  "mean: B - A" ~ b_GroupB) %>% mutate(model = 1),
  hdi(Post2, pars =  "mean: B - A" ~ b_GroupB) %>% mutate(model = 2),
  hdi(Post1, pars =  "mean: B - C" ~ (b_GroupB - b_GroupC)) %>% mutate(model = 1),
  hdi(Post2, pars =  "mean: B - C" ~ (b_GroupB - b_GroupC)) %>% mutate(model = 2)
) %>% select(model, contrast = par, lo, hi, prob)

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

##   model      contrast        lo         hi prob
## 1     1 "mean: B - A" -1.217038  5.1443755 0.95
## 2     2 "mean: B - A" -1.365880  5.1951917 0.95
## 3     1 "mean: B - C" -6.212982  0.2230876 0.95
## 4     2 "mean: B - C" -3.623713 -2.4066558 0.95

Comparing via LOO or WAIC shows that Model 2 is to be strongly preferred over Model 1. There is very little difference in the WAIC vs LOO estimates of elpd.

library(loo)

## This is loo version 2.1.0.9000.
## **NOTE: As of version 2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. Visit mc-stan.org/loo/news for details on other changes.

loo_compare(loo(model1_brm), loo(model2_brm))

##            elpd_diff se_diff
## model2_brm   0.0       0.0  
## model1_brm -64.5       8.9

loo_compare(waic(model1_brm), waic(model2_brm))

##            elpd_diff se_diff
## model2_brm   0.0       0.0  
## model1_brm -64.6       8.9

fert4_brm <-
  brm(Yield ~ Till * Fert + (1 | Field), data = SplitPlotAgri)

## Compiling the C++ model

## Start sampling

fert4a_brm <-
  brm(Yield ~ Till * Fert + (1 | Field), data = SplitPlotAgri,
      family = student())

## Compiling the C++ model
## Start sampling

Both WAIC and LOO find the two models to be very similar in estimated out-of-sample prediction accuracy. (The differences are about the same size as the standard errors, so no strong support for one model over the other.)

loo_compare(loo(fert4_brm, reloo = TRUE), loo(fert4a_brm, reloo = TRUE))

## 8 problematic observation(s) found.
## The model will be refit 8 times.

## 
## Fitting model 1 out of 8 (leaving out observation 15)

## 
## Fitting model 2 out of 8 (leaving out observation 20)

## 
## Fitting model 3 out of 8 (leaving out observation 31)

## 
## Fitting model 4 out of 8 (leaving out observation 73)

## 
## Fitting model 5 out of 8 (leaving out observation 77)

## 
## Fitting model 6 out of 8 (leaving out observation 78)

## 
## Fitting model 7 out of 8 (leaving out observation 83)

## 
## Fitting model 8 out of 8 (leaving out observation 97)

## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling
## Start sampling

## No problematic observations found. Returning the original 'loo' object.

##            elpd_diff se_diff
## fert4a_brm  0.0       0.0   
## fert4_brm  -1.4       0.9

loo_compare(waic(fert4_brm), waic(fert4a_brm))

##            elpd_diff se_diff
## fert4_brm   0.0       0.0   
## fert4a_brm -1.1       1.0